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This paper presents a coupled vortex-lattice flight dynamic model with an aeroelastic finite-element model 
to predict dynamic characteristics of a flexible wing transport aircraft. The aircraft model is based on NASA 
Generic Transport Model (GTM) with representative mass and stiffness properties to achieve a wing tip de- 
flection about twice that of a conventional transport aircraft (10% versus 5%). This flexible wing transport 
aircraft is referred to as an Elastically Shaped Aircraft Concept (ESAC) which is equipped with a Variable 
Camber Continuous Trailing Edge Flap (VCCTEF) system for active wing shaping control for drag reduction. 
A vortex-lattice aerodynamic model of the ESAC is developed and is coupled with an aeroelastic finite-element 
model via an automated geometry modeler. This coupled model is used to compute static and dynamic aeroe- 
lastic solutions. The deflection information from the finite-element model and the vortex-lattice model is used 
to compute unsteady contributions to the aerodynamic force and moment coefficients. A coupled aeroelastic- 
longitudinal flight dynamic model is developed by coupling the finite-element model with the rigid-body flight 
dynamic model of the GTM. 


I. Introduction 

The aircraft industry has been responding to the need for energy-efficient aircraft by redesigning airframes to be 
aerodynamically efficient, employing light-weight materials for aircraft structures, and incorporating more energy- 
efficient aircraft engines. Reducing airframe operational empty weight (OEW) using advanced composite materials 
is one of the major considerations for improving energy efficiency. Modern light-weight materials can provide less 
structural rigidity while maintaining sufficient load-carrying capacity. As structural flexibility increases, aeroelastic 
interactions with aerodynamic forces and moments can alter aircraft aerodynamics significantly, thereby potentially 
degrading aerodynamic efficiency. 

Under the Fundamental Aeronautics Program at the NASA Aeronautics Research Mission Directorate, the Fixed 
Wing project is conducting multidisciplinary foundational research to investigate advanced concepts and technologies 
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for future aircraft systems. A NASA study entitled “Elastically Shape Future Air Vehicle Concept” was conducted in 
2010 1 to examine new concepts that can enable active control of wing aeroelasticity to achieve drag reduction. This 
study showed that highly flexible wing aerodynamic surfaces can be elastically shaped in-flight by active control of 
wing twist and vertical deflection in order to optimize the local angle of attack of wing sections to improve aerodynamic 
efficiency through drag reduction during cruise and enhanced lift performance during take-off and landing. 

The study shows that active aeroelastic wing shaping control can have a potential drag reduction benefit. Conven- 
tional flap and slat devices inherently generate drag as they increase lift. The study shows that conventional flap and 
slat systems are not aerodynamic ally efficient for use in active aeroelastic wing shaping control for drag reduction. A 
new flap concept, referred to as Variable Camber Continuous Trailing Edge Flap (VCCTEF) system, was conceived 
by NASA to address this need. 1-2 Initial results indicate that the VCCTEF system may offer a potential pay-off for 
drag reduction that will result in significant fuel savings. In order to realize the potential benefit of drag reduction 
by active aeroelastic wing shaping control, configuration changes in high-lift devices have to be a part of the wing 
shaping control strategy. 

NASA and Boeing are currently conducting a joint study to develop the VCCTEF further under the research 
element Active Aeroelastic Shape Control (AASC) within the Fixed Wing project. 3-4 This study built upon the devel- 
opment of the VCCTEF system for NASA Generic Transport Model (GTM) which is essentially based on the B757 
airframe, 5 employing light-weight shaped memory alloy (SMA) technology for actuation and three separate chord- 
wise segments shaped to provide a variable camber to the flap. This cambered flap has potential for drag reduction as 
compared to a conventional straight, plain flap. The flap is also made up of individual 2-foot spanwise sections which 
enable different flap setting at each flap spanwise position. This results in the ability to control the wing twist shape 
as a function of span, resulting in a change to the wing twist to establish the best lift-to-drag ratio (L/D) at any aircraft 
gross weight or mission segment. Current wing twist on commercial transports is permanently set for one cruise con- 
figuration, usually for a 50% loading or mid-point on the gross weight schedule. The VCCTEF offers different wing 
twist settings for each gross weight condition and also different settings for climb, cruise and descent, a major factor 
in obtaining best L/D conditions. 

The second feature of VCCTEF is a continuous trailing edge flap. The individual 2-foot spanwise flap sections 
are connected with a flexible covering, so no breaks can occur in the flap platforms, thus reducing drag by eliminating 
these breaks in the flap continuity which otherwise would generate vorticity that results in a drag increase and also 
contributes to airframe noise. This continuous trailing edge flap design combined with the flap camber result in lower 
drag increase during flap deflections. In addition, it also offers a potential noise reduction benefit. 

The VCCTEF serves multiple functions as: 

• A wing shaping control device to twist the flexible wing to obtain changes in lift-to-drag ratios that will reduce 
cruise drag throughout the flight envelope. 

• A high-lift device for take-off, climb-out, let-down and final approach by using the full span cambered flap. 

• A full span roll control effector in lieu of traditional ailerons using the aft section of the cambered flap. 

• An aeroservoelastic (ASE) control device to compensate for reduced flutter margins of flexible wings. 

This paper describes an aeroelastic formulation of a flexible wing aircraft based on a one-dimensional structural 
dynamic theory that models the wing structure as a beam in a coupled bending-torsion motion. The aeroelastic angle 
of attack is derived from kinematics of aircraft rigid-body velocities and wing structural deflection velocities. The 
resulting nonlinear aeroelastic equations of bending-torsion motion are coupled with the aircraft rigid-body flight 
dynamic equations of motion. The nonlinear aeroelastic formulation takes into account the engine thrust forces which 
are coupled with wing aeroelasticity as a force follower. The formulation therefore is an aero-propulsive-elasticity. 
This inclusion of the propulsive effect may be important for aircraft with engine-mounted high flexible wing structures. 
The aeroelastic analysis also takes into account wing pre-twist and dihedral angles which can cause a high degree of 
coupling between the wing aeroelastic deflections and the aircraft rigid-body motion. 

In the present study, the aeroelasticity equations are transformed into a system of aeroelastic state-space equations 
using the finite-element method (FEM), a powerful numerical technique which converts a given system of partial dif- 
ferential equations into a truncated, discretized weak-form solution formulation which utilizes locally-defined basis 
functions (typically polynomials) to numerically approximate the solution to the governing partial differential equa- 
tions. In general, the standard finite-element method belongs to a class of numerical techniques known as weighted- 
residual methods, such as the Galerkin method, which seek to minimize the error between the “true” solution and 
the approximation space of basis functions. Mathematically, this property arises from the fact that the finite-element 
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method solves a variational weak form of the initial boundary value problem using arbitrary test functions which be- 
long to the Hilbert space of functions which are square-integrable. Proper application of the finite-element method 
will result in the creation of a matrix system of equations which may be solved using standard numerical techniques. 
A modal analysis is then performed to assess aeroelastic stability of the aircraft. Frequencies and damping ratios of 
the symmetric and anti-symmetric modes are computed. 

A vortex-lattice model of the flexible wing GTM, herein referred to as Elastically Shaped Aircraft Concept 
(ESAC), 6 is developed for coupling with the flight dynamic model and the aeroelastic finite-element model (FEM) 
to provide stability and control derivatives for the flight dynamic model and the aerodynamic loading information for 
the FEM. An automated geometry modeler developed in MATLAB provides an update of the aircraft wing deformed 
geometry for the vortex-lattice model. 

II. Description of Elastically Shaped Aircraft Concept 



Fig. 1 - Generic Transport Model (GTM) and Remotely Piloted Vehicle at NASA Langley 



Fig. 2 - GTM Planform 

The elastically shaped aircraft concept is modeled as a notional single-aisle, mid-size, 200-passenger aircraft. The 
geometry of the ESAC is obtained by scaling up the geometry of NASA generic transport model (GTM) by a scale of 
200: 1 1 . The GTM is a research platform that includes a wind tunnel model and a remotely piloted vehicle, as shown 
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in Fig. I. 4 5 Figure 2 is an illustration of the GTM planform. The reason for selecting the GTM is that there already 
exists an extensive wind tunnel aerodynamic database that could be used for validation in the study. The benchmark 
configuration represents one of the most common types of transport aircraft in the commercial aviation sector that 
provides short-to-medium range passenger carrying capacities. 

The aircraft has a take-off weight of 200,000 lbs for a typical operating load (gear up, flap up) that includes cargo, 
fuel, and passengers. Fuel weighs about 50,000 lbs for a range of about 3,000 nautical miles. 

To compute the mass and inertia properties of the benchmark aircraft, a component-based approach is used. The 
aircraft is divided into the following components: fuselage, wings, horizontal tails, vertical tail, engines, operational 
empty weight (OEW) equipment, and typical load including passengers, cargo, and fuel. The fuselage, wings, horizon- 
tal tails, and vertical tail are modeled as shell structures with constant wall thicknesses. 7 Based on publicly available 
data of component weight breakdown for various aircraft, 8 an average wing mass relative to the total empty weight of 
the aircraft is taken to be 24.2% of the OEW. 

To enable active wing shaping control, the wing structures of the ESAC are designed to increase wing flexibility. 
The wing bending and torsional stiffness quantities are designed to achieve a wing deflection that is about double 
of that of a conventional aircraft wing. The VCCTEF is divided into 14 sections attached to the outer wing and 3 
sections attached to the inner wing, as shown in Fig. 3. Each 24-inch section has three camber flap segments that can 
be individually commanded, as shown in Fig. 4. These camber flaps are joined to the next section by a flexible and 
supported material (shown in blue) installed with the same shape as the camber and thus providing continuous flaps 
throughout the wing span with no drag producing gaps. The flexible skin materials that cover the spanwise camber flap 
sections create constraints to the flap deflections. These constraints impose a certain relative flap deflection between 
any two adjacent spanwise flap sections. 

Using the camber positioning, a full-span, low-drag, high-lift configuration can be activated that has no drag 
producing gaps and a low flap noise signature. This is shown in Fig. 5. To further augment lift, a slotted flap 
configuration is formed by an air passage between the wing and the inner flap that serves to improve airflow over the 
flap and keep the flow attached. This air passage appears only when the flaps are extended in the high lift configuration. 

Because the wings are highly flexible, flutter margins can become a potential issue. Flight dynamics and control 
of a highly flexible wing aircraft must fully account for the effects of aeroelasticity. 



Fig. 3 - GTM with VCCTEF 


Variable Camber Flap 
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Fig. 5 - Cruise and High Lift VCCTEF Configurations 


III. Aeroelastic Analysis 

A. Reference Frames 



Fig. 6 - Aircraft Reference Frames 

Figure 6 illustrates three orthogonal views of a typical aircraft. Several reference frames are introduced to facilitate 
the rigid-body dynamic and structural dynamic analysis of the lifting surfaces. For example, the aircraft inertial 
reference frame A is defined by unit vectors ai, a2, and a 3 fixed to the non-rotating earth. The aircraft body-fixed 
reference frame B is defined by unit vectors b i , bo, and b3 aligned with the roll, pitch, and yaw axes, respectively. The 
right wing elastic reference frame C is defined by unit vectors Ci , C2, and C3. The reference frames B and C are related 
by three successive rotations: 1) the first rotation about b 3 by the sweep angle | + A of the elastic axis that results in 
an intermediate reference frame B defined by unit vectors b,, b,, and b 3 (not shown), 2) the second rotation about 

n / 

b 0 by the dihedral angle T of the elastic axis that results in an intermediate reference frame C defined by unit vectors 
Cj, C9- and c 3 (not shown), and 3) the third rotation about c, by an angle j r that results in the reference frame C. This 
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relationship is expressed as 
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The left wing elastic reference frame D is defined by unit vectors di , do, and d 3 . The reference frames B and D are 
related by three successive rotations: 1) the first rotation about b 3 by the elastic axis sweep angle f + A that results 
in an intermediate reference frame B ' defined by unit vectors bj, b^, and b 3 (not shown), 2 ) the second rotation about 
b" by the elastic axis dihedral angle F that results in an intermediate reference frame D defined by unit vectors dj , 
di, and d 3 (not shown), and 3) the third rotation about d[ by an angle n that results in the reference frame D. This 
relationship is expressed as 
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The aircraft velocity at the aircraft center of gravity (CG) in the aircraft body-fixed reference B is expressed in the 
reference frames C and D as 
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where (u.v,w) are the aircraft velocity components in the forward, lateral, and downward directions defined by the 
unit vectors (b 1 , bi , b 3 ), respectively. 

Generally, the effect of the dihedral angle can be significant. In the analysis, the aeroelastic effects on the fuselage, 
horizontal stabilizers, and vertical stabilizer are not considered, but the analytical method can be formulated for ana- 
lyzing these lifting surfaces if necessary. In general, a whole aircraft analysis approach should be conducted to provide 
a comprehensive assessment of the effect of structural flexibility on aircraft performance and stability. However, the 
scope of this study pertains only to the wing structures. 


B. Elastic Analysis 

In the subsequent analysis, the combined motion of the left wing is considered. The motion of the right wing is a 
mirror image of that of the left wing for symmetric flight. The wing has a varying pre-twist angle y(x) commonly 
designed in many aircraft. Typically, the wing pre -twist angle varies from being nose-up at the wing root to nose-down 
at the wing tip. The nose-down pre-twist at the wing tip is designed to delay stall onsets. This is called a wash-out 
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twist distribution. Under aerodynamic forces and moments, the aeroelastic deflections of a wing introduce stresses 
and strains into the wing structure. The internal structure of a wing typically comprises a complex arrangement of load 
carrying spars and wing boxes. Nonetheless, the elastic behavior of a wing can be captured by the use of equivalent 
stiffness properties. These properties can be derived from structural certification testing that yields information about 
wing deflections as a function of loading. For high aspect ratio wings, an equivalent beam approach can be used to 
analyze aeroelastic deflections with good accuracy. The equivalent beam approach is a typical formulation in many 
aeroelasticity studies. 9 It is assumed that the effect of wing curvature is ignored and the straight beam theory is used 
to model the wing deflection. The axial or extensional deflection of a wing is generally very small and therefore can 
usually be neglected. 

Consider an airfoil section on the left wing as shown in Fig. 7 undergoing bending and torsional deflections. Let 
(x,y,z) be the coordinates of point Q on a wing airfoil section. Then the undeformed local airfoil coordinates of point 
Q are 


y 
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sin y cos y 
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where T) and t, are local airfoil coordinates, and y is the wing section pre -twist angle, positive nose-down. 10 
Then differentiating with respect to x gives 
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Fig. 7 - Left Wing Reference Frame of Wing in Combined Bending-Torsion 


Let © be a torsional twist angle about the A-axis, positive nose-down, and let W and V be flapwise and chordwise 
bending deflections of point Q, respectively. Then, the rotation angle due to the elastic deformation can be expressed 
as 

0 (x,t) = ©di - W v d 2 + V v d 3 (7) 

where the subscripts x and t denote the partial derivatives of 0, W, and V. 

Let (xi,yi,zi) be the coordinates of point Q on the airfoil in the reference frame D and p be its position vector. 
Then the coordinates (xi,yi,zi) are computed using the small angle approximation as 
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Neglecting the transverse shear effect, the longitudinal strain is computed as 1 
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where 
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where E is the Young’s modulus; G is the shear modulus; y is the derivative of the wing pre -twist angle; I yy , I yz , and 
/,, are the section area moments of inertia about the flapwise axis; J is the torsional constant; and B \ , B 2 , and B 3 are 
the bending-torsion coupling constants which are defined as 
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The strain analysis shows that, for a pre-twisted wing, the bending deflections are coupled to the torsional deflection 
via the slope of the wing pre-twist angle. This coupling can be significant if the wash-out slope y is dominant for 
highly twisted wings. 


C. Aeroelastic Angle of Attack 

The relative velocity of the air approaching a wing section includes the contribution from the wing elastic deflection 
that results in changes in the local angle of attack. Since aerodynamic forces and moments are dependent on the local 
angle of attack, the wing aeroelastic deflections will generate additional elastic forces and moments. The local angle 
of attack depends on the relative approaching air velocity as well as the rotation angle 0 from Eq. (7). The relative 
air velocity in turn also depends on the deflection-induced velocity. The velocity at point Q due to the aircraft velocity 
and angular velocity in the reference frame D is computed as 

v e =v + fflxr = (ubi + vb 2 + wb 3 ) + (pb t + qb 2 + rb 3 ) x (-x fl bi - y„b 2 - z„b 3 ) 

= (u + ry a - qza) bi + (v - rx a + pz a ) b 2 + (w + qx a - py a ) b 3 = x t A 1 + y, d 2 + z t d 3 (18) 
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and ( p,q,r ) are aircraft angular velocity components in the roll, pitch, and yaw axes, and (x a , y a , z a ) is the coordinate 
of point Q in the aircraft body-fixed reference frame B relative to the aircraft C.G. (center of gravity) such that x a is 
positive when point Q is aft of the aircraft CG, y a is positive when point Q is toward the left wing from the aircraft 
C.G., and z a is positive when point Q is above the aircraft C.G. 

The local velocity at point Q due to aircraft rigid-body dynamics and aeroelastic deflections in the reference frame 
D is obtained as 10 
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In order to compute the aeroelastic forces and moments, the velocity must be transformed from the reference frame 
D to the airfoil local coordinate reference frame defined by (p, rj, |) as shown in Fig. 2. Then the transformation can 
be performed using successive rotation matrix multiplication operations as 
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for small deflections. 

The local aeroelastic angle of attack on the airfoil section due to the velocity components Vjj and v<e in the reference 
frame D, as shown in Fig. 8, is computed as 
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where w, is the local downwash due to three-dimensional lift distribution over a finite wing and 
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Assuming that the wing dihedral T is small and neglecting the local downwash w,, the local aeroelastic angle of 
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attack is evaluated as 
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Fig. 8 - Airfoil Local Coordinates 


The aeroelastic angle of attack is generally a nonlinear function of the structural deflections. A nonlinear aeroelas- 
tic theory based on this approach has been developed. 12 


In this study, a linear aeroelasticity is developed by neglecting nonlinear deflection-dependent terms and letting 
u ~ Voc, v « vjp, and w ~ VLa. Thus 


(Voo + ry a -qz a ) sin AT + ( VU J3 - rx a + pz a ) cos Ar - ( V«, a + qx a - py a ) + W, + y 0, 


+ 


V„o cos A 

[~ ( y °° + r 7a ~ <lZa) sin A ~ (V°°P - rx a + PZg) COS A - (Too 0! + qx a - py a ) T] W x 

VLcosA 


+ 


— (Voo + ry a - qz a ) cos A + (Voo/3 - rx a + pz a ) sin A] (0 + 7) 

(A, cos A 

(Voo + ry a - qzg ) sin AT + (Vooj3 - rx a + pz a ) cos Ar - (TopOl + qx a - py a ) 

V 2 cos 2 A 


x 


x 


{rya - qz a ) COS A + (Voo/3 - rx a + pz a ) sin A + V, - z& t 


+ [(Voo + ry a - qza) sin A + (Vooj3 - rx a + pz a ) cos A + (V»a + qx a - py a ) T] V x 


+ [(Voo + ry a - qza) sinAr + (V»J3 - rx a + pz a ) cos Ar - (Voo a +qx a - py a )\ (© + 7) 


(28) 
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Assuming z ~ 0, then evaluating the various partial derivatives of the local aeroelastic angle of attack yields 
(Xq = — y — tanAr(l + tan Ary) 
da c l+2tanAry 
da cosA 

=tanAy— r(sec 2 A + 2tanAry) 
dp 

da c y a (1 +2tanAry) z a (sec 2 Ar~tanAy+2tanAr 2 y) 

dp VcoCosA Voc 

da c _x a (l+2tanAry) z a J (l +2tan 2 Ar 2 ) 
d q Voc cos A Voc 

da c x a (rsec 2 A — tanAy+2tanAr 2 y) y a y (l +2tan 2 Ar 2 ) 

~d7 = v: v: 

da, _ y 

da 2 cos 2 A 

da c ^ ^ , 

=-r(tanA + ry) 

da, _ yjy z 2 r(tanA + Ty) y a z a (tanA + 2ry) 

dp 2 V^cos 2 A V 2 Vi cos A 

da, x 2 y z 2 tanAr(l — tanAry) x a z a (1 — 2tanAry) 

dq 2 Vicos 2 A+ Vi Vicos A 

da c x 2 r(tanA + Ty) y 2 tanAr(l -tanAry) x a y a r (l -tan 2 A-2tanAry) 

IfA ~ yi 1 yi y2 

da c tanA + 2ry 
dafi cosA 

da c _ 2 y a y z a (tanA + 2ry) 
dap Vx, cos 2 A Voc cos A 

da, _ 2 x a y | z a (1 — 2tanAry) 

<9ag Voc cos 2 A Voc cos A 

da c x a (tanA + 2ry) y a (1 — 2tanAry) 

dar VooCosA V*,cosA 

da c y a (tanA + 2ry) 2z a r(tanA + Ty) 

dfip VooCosA Voc 

da c x a (tanA + 2ry) z a r(l -tan 2 A-2tanAry) 

dpq VocCOsA Vx, 

da c _ 2x a r (tan A + Ty) | y fl r(l -tan 2 A-2tanAry) 

w>- = v: + ^ 

da c z 2 r(l -tan 2 A-2tanAry) | 2 x a y a y | x Q z a (tanA + 2ry) y a z a (1 -2tanAry) 

dpq Vi Vicos 2 A Vicos A Vicos A 


T a 


(1 —2 tanAry) x a y a (tanA + 2ry) 2x a z a r (tan A + Yy) y a z a r (l -tan 2 A-2tanAry) 


dpr Vicos A Vicos A V 2 V 2 

da c _ x 2 (tanA + 2ry) x a y a (1 -2 tanAry) x fl z a r(l -tan 2 A-2tanAry) 
dgr VicosA VicosA V 2 

2y a Za tan AT (1 — tanAry) 

vi 

da, _ (Voo + ry a - qzg) sin A + (Voo/3 - rx a + pz a ) cos A + (VocOi + qx a - py a )T 
dW x ~ VicosA X 

", , (Voo + rya-qz a )smAry+(Voop-rx a +pZa)cosAry-(Vo 0 a + qx a -pya)y 

VccCosA 
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( 51 ) 


da c __ 1 _ (Voc + ry a - qzg) sin Ary + (Kofi - rx a + pz a ) cos Ary - (V^a+qxg - py a ) y 

dW, Voc cos A V^cos 2 A 

da c _ (Voc + ry a - qzg) sin A + (Vx,/3 - rx a + pz a ) cos A + (Voc a + qx a - py a ) F 

dV x Voc cos A 

T _ (Voc + ry a - qzg) sinAr + (V*,/) - rx a + pz a ) cos Ar - (VcO! + qx a - py a ) ~ 

X y VocCosA 

da c _ y _ (Voc + ry a - qzg ) sinAr + (V43 - rx a + pz a ) cos Ar - (V^a + qx a - py a ) 

<9V VocCosA V ( 2 cos 2 A 

^a c - _ (Voc + ry a - qza) cos A - (Voc/3 - + pz a ) sinA 

<9© VocCosA 

_ [(Voc + ry a - qzg) sinAr + (Voc)3 - rx a + pz a ) cos Ar - (Voca + qx a - py a )] 2 

Vi cos 2 A 

da c _ _ y _ yy [(Voc + ry a - qz a ) sinAr + (Voo/3 - rx a + pz a ) cos Ar - (Voca + qx a - py a )\ 
d& t VocCosA V ( 2 cos 2 A 


(52) 

(53) 


(54) 

(55) 


These partial derivatives contribute to the aeroelastic angle of attack as follows: 


, x da c 

0Cc{x,y,z) = ao+-^-s- 


da c . 

dW 3 


■W t 


da r. 


; W, 


where s = 


a 


P q r 


a “ 


q 


dW t 

~ a/3 ap aq 


<9a c . da CTr 

M Vx+ M Vt ' 


da c 
d © 


© 


da c 

d& t 


vector of the aircraft flight dynamic state variables. 

The terms W r , V x , and © contribute the aerodynamic stiffness, while the terms W t , V t , 
aerodynamic damping. 

The local angle of attack of an airfoil section at the aerodynamic center is evaluated at x a 
y = —e, and z = 0, where x ac is the forward distance of the aircraft center of gravity from the aerodynamic center of a 
wing section and e is the forward distance of the aerodynamic center from the elastic axis. Then 


■©, 


(56) 



T 

pq 

pr qr 

is a 

& t 

contribute 

to the 

ac-> ya — yac-> Za 

: Zac j 


a ar — ao T 


da rh 


da r „- , 


ds 


dU \ 


U x 


dOa, 

dU t 


U t - 


d a„ r d a ac „ 


<9© 


d&t 


(57) 


There is another source of lift, called non-circulatory lift due to the reaction force of the air volume surrounding a 
wing section. The non-circulatory lift force is based on the aeroelastic angle of attack at the mid-chord location which 
is evaluated at x a = x mc , y a = y mc , z a = Z m c, y = e m , and z — 0, where x mc is the forward distance of the aircraft center 
of gravity from the mid-chord location of a wing section and e m is the forward distance of the aeroelastic center from 
the mid-chord location. Then 


d a mc d a, nc 

Omc — Oq H ^ s H . , U x 


ds 


dUr 


d a„„ 
~dU t 


d a mr d a,;,,- 

U t + ^U-&+^L&, 


d& 


d&. 


(58) 


D. Aeroelastic Equations of Coupled Bending-Torsion Motion 

The equilibrium conditions for bending and torsion are expressed as 11 

(59) 

(60) 
(61) 

where m x is the pitching moment per unit span about the elastic axis, f z and f y are the lift and drag forces per unit 
span, respectively, and m y and m z are the bending moments per unit span about the flapwise and chordwise axes of the 
wing. 

The wing section lift coefficient is given by 

CL ac = c La C (k) a ac + c Ls S (62) 


dM x 

= - m x 
dx 

d 2 My duly 


<9 2 M, 

dx 2 


_ dm, 

dx 
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where k = is the reduced frequency parameter, a> is the frequency of wing oscillations, c is the section chord, 
cl u is the section lift curve slope, cl 5 is a vector of the section lift derivatives due to the VCCTEF deflection 5 = 

r it 


<5i <?2 ■ ■ • <5i4 

The function C(k ) is the Theodorsen’s complex-valued function 13 which is also expressed in terms of Bessel 
functions as 


C{k) =F(k) + iG(k) 


(63) 


where F (k) > 0 and G (k) < 0 are shown in Fig. 9. 

When k = 0, the airfoil motion is steady and C (k) is real and unity. 



Fig. 9 -Theodorsen’s Function 

The wing section lift coefficient due to harmonic motions is expressed as 

c Lac = c La a ac F (k) + c La + c Ls 8 (64) 

In addition, the apparent mass of the air contributes to the lift force acting at the mid-chord location as follows: 

KCC mc C 




2Voc 


(65) 


The total section lift coefficient is 


c l — cl„ c + c Lmc 


The section pitching moment coefficient is evaluated as 


. e e m . 

Cm — C mac + c CL ac ^ CL mc + C mg O 


(66) 


(67) 


where c mac is the section pitching moment coefficient at the aerodynamic center and c mg is a vector of the section 
pitching moment derivatives at the elastic axis due to the VCCTEF deflection. 

The section drag coefficient is expressed in a parabolic drag polar form as 


cd = cd 0 + 


kARe 


(68) 


where cd 0 is the section parasitic drag coefficient, AR is the wing aspect ratio, and £ is the span efficiency factor. 
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Expanding the expression and ignoring the nonlinear terms for the aeroelastic analysis, the section drag coefficient 
is obtained as 


CD = CD o + CD a ClacF (k) + Co a 


a ac c G (k) 

2V.O k 


+ c D s 8 


(69) 


where 

2 c L a ° 

^ La OCo 

C ° s n A Re 


(70) 

(71) 


In addition, the propulsive effects of the aircraft engines must be accounted for in the analysis. Both the engine 
mass and thrust can contribute to the aeroelasticity . 12 The propulsive force and moment vector are computed as 


i e =8{x-x e ) 


T 0 m e g 


— sinA —cos A sinAr 


" d! ' 


(-T sinA — m e gT) di 

— cos A sinA cosAr 


d 2 

II 

Oo 

1 

X 

— T cos Ad 2 

1 

n 

0 

1 


^3 


(T sinAr — m e g ) d 3 


(72) 


m e =r e x f e = (x e di -y e &2 -z e d 3 ) x i e = 8 (x-x e ) 


(— Ty e sin AT — Tz e cos A + m e gy e ) di 
(—Tx e sin Ar + m e gx e + Tz e sinA + m e gz e r) di 
(— 7X,cosA — Ty e sinA — m e gy e T)d^ 


(73) 


where T is the engine thrust, m e is the engine mass, ( pc e ,y e ,z e ) is the coordinate of the engine thrust center such that y e 
is positive forward of the elastic axis and z e is positive below the elastic axis, and 8 (x — x e ) is the Dirac delta function 
such that 

j 8(x-x e )f(x)dx = f(x e ) (74) 

Transforming into the local coordinate reference frame and neglecting nonlinear contributions, the propulsive 
forces and moments are given by 


fx e = 8 (x-x e ) [-T sinA -m e gT-T cos AV X + {T sinAr- m e g)W x ] (75) 

fy e = 8 (x-x e ) [(T sin A + m e gT) V x - T cos A + (7 sinAr- m e g) (© + y)] (76) 

f Ze = 8 (x — x e ) [(T sin A + m e gr) W x + T cosA(0 + y) + T sin Ar — m e g\ (77) 

m Xe = 8 (x — x e ) [— Ty e sinAr — Tz e cos A + m e gy e + (—Tx e sinAr + m e gx e + Tz e sinA + m e gz e r) V x 

— (Tx e cosA+ Ty e sinA + m e gy e T) W x ] (78) 

>riy e = 8 (x — x e ) [(Ty e sinAr + Tz e cosA — m e gy e )V x — Tx e sin Ar + m e gx e + Tz e sin A + m e gz e F 

- (Tx e cos A + Ty e sinA + m e gy e T) (0 + 7 )] (79) 

m Ze = 8 (x — x e ) [( Ty e sinAr + Tz e cos A — m e gy e )W x — (—Tx e sinAr + m e gx e + Tz e sinA + m e gz e r) (0 + 7 ) 

—Tx e cosA— Ty e sin A — m e gy e T\ (80) 


The partial derivatives of the moment components are 


dm e x 

r)x 

dmi 

y_ 

dx 

dm l 

dx 


8 (x — x e ) 
8(x-x e ) 
8 (x-x e ) 


\{-Tx e sinAr + m e gx e + Tz e sinA + m e gz e H V xx - (T x e cos A + Ty e sin A + m e gy e T) W xx ] 

( Ty e sinAr + Tz e cos A - m e gy e ) V xx - (Tx e cos A + Ty e sinA + m e gy e T) (& x + 7 '^ 
\(Ty e smAr + Tze cos A - m e gy e ) W xx + (Tx e sin Ar - m e gx e -Tz e sin A - m e gz e T ) ( 0, + 7 


(81) 

(82) 


(83) 


Using the sign convention as shown in Fig. 10 the lift and drag forces and pitching moment per unit span can be 
expressed as 
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Fig. 10 - Airfoil Forces and Moment 


. . . c G(k) nce m . 

cc m ac + CCl 0 +eCL a (%acF \k) ~\~ €CL a CCac T , T ^ me 


qoo cos 2 Ac + mge cg — mk 2 ® tt 


fy = 


fz = 


2Voo k 2Voo 

+ me cg W tt - 8 (x - x e ) \m e (; y 2 + z 2 e ) ®, t - m e y e W„ + m e z e V t , ] + 8(x-x e ) [~Ty e sin Ar - Tz e cos A + m e gy e 
+ (— 7A e sinAr + m e gx e + 7 T z e sinA-)-m f ,gz ( >r) V x — (Tx e cosA. + Ty e sinA + m e gy e T)W x \ (84) 

c G{k)' 


C D 0 + CDi + CD a a acF (k) + Cj) a a ch 


TVoo k 

+ S (x — x e ) [(r sinA + w e gr) V* — T cosA+ (r sin Ar — m e g) (0 + 7)] 


qoo cos 2 Ac — mV tt + 8 (x — x e ) {—m e V tt — m e z e ®tt) 


(85) 


. c G(k) nc . 

CLq + CL a CC ac F (k) + CL a CC ac — - h (X mc 


qoo cos' Ac — mg — mW n + me cg @tr 


+ 8 (x — x e ) (— m e W t t +m e y e ® tt ) + 8 (x — x e ) [(T sinA + m e gr)lF v + r cosA(0 + 7 ) +T sinAr — m e g] ( 86 ) 

where q 0 0 is the dynamic pressure, m is the wing mass distribution, e cg is the eccentricity between the center of mass 
and the elastic axis (positive corresponding to the center of mass located forward of the elastic axis), k is the torsional 
radius of gyration, and the term cos 2 A accounts for the wing sweep angle A as measured from the elastic axis. 

The bending and torsion aeroelastic equations then become 

d 

dx 


GJ + EB X (y) 
cc >n ac + ec La a ac F (k) + ec La a ac 


® x -EB 2 yW xx -EB 3 yV xx = 
c G(k) nce m . 


2Voo k 2Voo 

,2 1 „2\ 


a,nc + (eCL s +cc,n s ) 8 


qoo cos 2 Ac — mge cg + mk 2 @ tt 


- me cg W„ + S(x-x e ) m e (y 2 + z 2 ) ©„ - m e y e W tt + m e z e V n + Ty e sinAr + 7> e cosA — m e gy e 

— (—Tx e sinAr + m e gx e + Tz e sin A + m e gz e F) V x + (7x e cos A + Ty e sin A + m e gy e T) W x 


(87) 


d 2 
dx 2 


( -EB 2 y® X + EI yy W XX - EIy Z V XX ) = 


. c G(k) nc . 

CL a &acF \k) + CL a CC ac ~ b 2V~ ° ^ 


qoo cos“ Ac — mg — mW„ + me cg ®tt 


+ 8 (x — x e ) — m e W tt + m e y e ® tt + ( T sin A + m e gT) W x + T cos A (0 + 7) + T sinAr — m e g 
- (Ty e sin Ar + T z e cos A - m e gy e ) V xx + (Tx e cos A + Ty e sin A + m e gy e T) ( © v + 7 


( 88 ) 
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d 2 
dx 2 


(-EB 3 y 


©x — EIy Z W XX + EI ZZ V XX ) — 


(: D„ T C D„_ (%acF (k) ~T C j) a a, l( 


c G(k) 


-c Dl .8 


q oo cos“ Ac — mV tt 


2Voo k ' ~' JS 

+ 8 {x — Xe) | —m e V tt — m e z e ®tt + (T sin A + m e gY) V x — T cos A + (T sin AT — m e g) (© + y) 

- (Ty e sin Ar + Tz e cos A - m e gy e ) W xx - (T x e sin Ar - m e gx e - Tz e sin A - m e gz e T) [© A - + yj 


( 89 ) 


Taking advantage of symmetry, the motion could be decomposed into symmetric and anti-symmetric motions 
subject to either symmetric-mode boundary conditions W x (O. t) = V x (O.t) = 0 or anti-symmetric mode boundary 
conditions 0(O,f) = W (0 ,t) = V (0 ,t) = 0 at the left end. Half of the mass and mass inertia of the aircraft structure 
without the wings are added to the generalized mass of the system. Defining the vector quantities 


U = 


co = 


Ca — 


eg = 


W 

V 

0 

CDo 

C La 

C D a 

C L S 

C D S 

TIC 

2Voo 

c 

8 
0 

e cg 


(90) 

(91) 

(92) 

(93) 

(94) 

(95) 

(96) 


r e 


fo 


fe 


f&x 


fu x 


ye 

~Ze 

T sinAr — m e g 
—T cosA 

T cos A 

T sin AT — m,,g 

Tx e cosA + Ty e sin A + m e gy e V 

— ( Tx e sinAr — m e gx e — Tz e sinA — m e gz e r) 

T x e cos A + Ty e sin A + m e gy e T 

— (—Tx e sinAr + m e gx e + Tz e ^A. + m e gz e ^') 


(97) 

(98) 

(99) 
(100) 
(101) 


B = 
1 = 

J = 


b 2 b 3 j 

Iyy —lyz 

~k z y h.z 

0 1 
1 0 


( 102 ) 

(103) 

(104) 
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The aeroelastic partial differential equations are given as 


|^| GJ + EBx^y ® x -EByU xx j = 

_ . C G (k') . / No 2 a 

'mac ^L 0 “1“ ^acF \k) €C L a CCac ^ ^ me H - H” CC m ^ ) C C[oo COS Ac 

- mge cg + mk 2 ® tt - me cg U,t + 8 (x - x e ) m e (y 2 + z 2 ) ® tt - m e rj U tt 

+Ty e sin Ar + Tz e cos A — m e gy e + fu x U x \ (105) 

d 2 ( T ' ^ \ 

^ y—EB y 0 V + EIU xx j = 

c G (k') 

co + c a CCacF ( k) + c a CCac yyy + c c a mc + c$8 qoo cos 2 Ac - ma 

- mUtt + me cg ®„ + 8(x-x e ) \-m e U n + m e r e ®„ + (T sin A + m e gT) U x 

+/o - (Ty e sin Ar + Tz,, cos A - m e gy e ) JU XX + /©(© + y) + f@ x (® x + 7 ) (106) 

IV. Finite Element Modeling 

A. Finite-Element Discretization 

The aeroelastic equations describe the wing bending and torsional deflections due to aerodynamic forces and moments. 
Using the finite-element method, 14 the structure can be discretized into n equally spaced one-dimensional elements. 
Then the bending and torsional deflections can be approximated as 

0(m) = ^0, (M) 

i= 1 

U (x,t)=iui(x,t) 

i= 1 

where i is the ;-th element, and n is the number of nodes. 

For each element, the bending and torsional deflections are approximated as 

0; (-M) = [ m (x) y /2 (x) | d }‘ y = Ne (x) 9j (t) 

L J y 2 , (t ) 


Wl i (0 

w 'lf (0 

Vli (0 

4, (0 

W2i (t) 

W 2i (0 
M (0 

. 4,.(0 

where the subscripts 1 and 2 denote values at nodes 1 and 2 , and xffj (x), j = 1 , 2 , and 0 * (x), k - 1 , 2 , 3, 4 are the linear 
and Hermite polynomial shape functions 

VO (x) = 1 — y ( m ) 

Vfc(*)=y (112) 
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v .( xt \ = ('iW fcW 0 0 04 0) o o 

0 0 01 (x) 02 (x) 0 0 03 (x) 04 (x) 


(109) 


= N u (x)ui(t ) (110) 


(107) 

(108) 
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<M*)=l-3(y) 2 +2(y) 3 (113) 

^( y ) 2 - 2 © 3 

«. 4 M=([-()) 2 + 0) 5 ] (116) 

where x £ [0,/] is the local coordinate and / = is the element length. 

It can be shown that Hermite cubic shape functions result in exact nodal displacements and slopes, making them 
ideal candidates for problems involving beam bending. 

The weak-form integral expressions of the dynamic aeroelastic equations are obtained by multiplying the bending 
and torsion aeroelastic equations by Nj (x) and Nj (x), and then integrating over the wing span. This yields 

£/a^{ GJ + EBi^y N g 0j — EBy N^uX dx = 
i=l 0 K l J ) 

V" 1 / a tT f ^ ^ ac d CC ac t d CC ac d 0C ac ATC . d 0C ac • \ . . 

L J N e \cc mac + ec La (a< ) + —s+ ~^N u Ui + ~^N uUi + — Ng 6, + —Ng 0.) F (k) 

f d CC a c . d (X a c ,/ . dct a c AT .. d 0C a c AT a dct a c AT C G (Ji) 

+eCLa Ku, ^ uUi ^ e * ^ e i ) 2v:~ 

f dctmc . d CCmc ^ CC m c .. dCL m c AT a dCL m c ( \ c 2 k j 

~ 214, \^r s + su x Nuli ‘ + -Ju; Ni<Ui + -ji & Nedi + d@, Ne0 j + ( eCL * +CCm ^ ^ cos Acdx 
« ' 

+ / Ng (—mge C g+mk 2 Nedj — m£ C gN u Ui)dx 

1=1 o 

n r -| 

+ T J Nj m e (y 2 e + z 2 ) N e 0i - m e r] N u iii + Ty e smAT + Tz e cosA - m e gy e + f Ux N' u Ui (117) 

i=i L 2x ~ Xe 

t ( K -EB T yN e e i +EIN" l u i ) dx = 

1=1 0 

/ 

f atT i | d(X ac doCac' d(X ac AT . . dcC ac AT a d(X a c A1 a\ t?(i\ 

5 J N- co + c a [ao+^rS+ N u u t + ^ N u u t + —N e 0 , + -^-N e 0, J F ( k ) 

1-1 o 

( da ac . dcc ac i . doc ac .. doc ac a da ac a \ c G (k) 

+c « (iiT + 7 ^ a, '“' + 7wr w '“ i+ w w » 9i+ <j 07 "» e 'J 

+c c (%*+ + + + +..'5«] 9-cor Ac* 

/ 

n r n r f 

+ ^ / ivj (—ma — mN„iii + m£ cg Ng 0,) Jx + ^ —m e N u Ui + m e r e Ng 9j + (T sin A + m e gr) A„m/ 
i=1 o ,=1 

+/o - (7>e sin Ar + rz c cos A - m e gy e ) Tiv" «,■ + /© (A^ 0 0; + 7) + /© t 0/ + 7 ) (118) 

The expressions of the left hand sides can be integrated by parts upon enforcing the boundary conditions as 

j N °^L { G/ + £B i(l ')" N e e i - EByNlu^dx = - jN 0 J j GJ + EBi (y') 2 A e 0,- -EByN'u^ dx (119) 
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i 2 1 

I Nj ^ ( K -EB T yN e e i +EIN",u^j dx = J iv" T (-EB T 'yN' e 6 i +EIN''u^j dx (120) 


o o 

Upon substitution, the final form of the aeroelastic equations are produced as 


M 

G/ + £Bi (y ) 2 

o 1 

L J 


L N e 

' =1 o 


N e ft — EBy N u Ui > dx 


d OCqc d OC a c 


da. 


d a nr , T ^ d a nr , 


+ ecL a ( Oq + ~^r s + N u Ui + -^-N u Ui + A^ft + - 37 ^-Afeft ) F (fc) 


<9F, 


. d a ac . d a ac f . da ac .. da ac • da ac 
+ eCLA — s+ — -N u ui + 0i + Ng 9j 


d® 

ftce m ( da mc . d a mc r da mc . da mc • da mc ~ \ / \ o 

^ — s + w w “"' + ^r N ""' + e ‘ + w " 8 e ‘ ) + ( " Li + > 5 


5© 01 5©, 

c G(£) 

2 V^~k~ 


qoo cos 2 Acc/x 


/ 

n /• 

+ (—mge cg +mk 2 NgQj — m£ cg N u iii)dx 

i=ii 


+ £ Nj Ule (^e + - 'W«rj Al„Mj + Ty e sin Ar + Tz e cos A — m e gy e + fu x N u Ui 
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The mass matrix, damping matrix, stiffness matrix, and force vector corresponding to each finite element i are then 
established as 


in 

N J e k 2 N e 

Nq £ C gN u 

dx + m e 

’ {y 2 e +£)Ne 

-NjrjN, , ' 

J 

0 

_ ~Nje cs N e 

NjN u 


- Njr e N e 



+ 


At I ( doCgc c G( A) Tlc'iti dcc, nc \ a ; T ( dcc ac c uc(' nl otx mc \ 

tve 5@f 2Ko j. 2V» 5®, y ■' V0 ■' v e \^ ec ic< ay, 2K» A 2V«, 5y, I JV « 


t T gofac C G(A) rtce m da, 

"e y cc L<x d®, 2Voo k 2V» 56 

AtT / ^Cfac c G(A) 5a mc A w W T / 5ofa c c G(fc) 

yv » ^ c « 5©7 2 i^^ +Cc_ 5©rj Ve v » v c “5gT2u:^ 


c G(A) tree,,, da m 


i r ^Otmc | iu 
+ 5£/f |tv„ 


quods' Acdx (123) 


19 of 38 


American Institute of Aeronautics and Astronautics 



Q = 


Nj ec La (^F(k) + ^^^]N e 1V 0 T ec ia ffiF (k) N u 


da, 


G(k )' 


ddcu 


- Njc a 


aeT^W- 


dttac c G(k) 
d® 2VU 1 


Ne -Njc a ^F(k)N„ 


qoo cos 2 A cdx 


0 N J eci d<Xac — G ^ n' 

V JV 0 ec L a dLJ x 2Voo k JV « 

0 — N T r - a “ c c N 

u yV i i c a d Ux 2V«, k iV « 


L 

q 00 cos 2 Acdx + J 


AT ~ r d(X mc AT „ AT I ACe m QLX mc at' 

iy e 2Voo de iy ° iy e 2Ko du x iV * 


T Ttce m doc mc at ' 


-Njc c ^N e 


— r — a '" c N 
Jv n Cc dU x ly u 


qoc cos" A cdx 
(124) 


K; 


N a 


GJ + EBi (j)‘ 

~n'' t eb t yn' 


N g -N e J EB T N u 

n" t ein" 


dx 


0 

NjfuX 


Nj(f®N e +f® x N' d ) 

(T sin A + m e gT)N' u ~{Ty e sin Ar + Tz e cos A - m e gy e ) J/v'' 



N J ec L a d -^F(k)N 9 Nj ec La d -^F ( k ) n' u 
-N j Ca^F (k) Ne —Nj c a ^F (k) N u 


cos 2 A cdx 


(125) 


Fi = m 


N eS e c g 
—Nj a 


dx + 


-Nj (: Ty e sin Ar + Tz e cos A - m e gy e ) 


dFi 

ds 

dFi 

ds 

dFi 

d8 


-Nj (cc mac + ec La a o) 
Nj (co + C a CCo) 

-N J e ec La d -^F(k) 
Njc La ^F(k) 


Nj (fo+feY+faY 


q^ cos 2 A cdx 


qoo cos 2 A cdx 


— N~I f PCT OUac c yJyK) — 

ly 9 \ ec La ds 2Voo k 


d cCqc c nce m dcc m 


2Voo ds 


xrT ( r dOCqc C G{k) | ,, d(Xmc 
iy u ds 2Voo k ~^ Cc ds 


qoo cos 2 A cdx 


l r 


0 L 


~ N e ( ec ^s +cc ms ) 


Njc s 


qoo cos 2 A cdx 


(126) 


(127) 


(128) 
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The discrete global approximation system is formed by enforcing equilibrium conditions at element interfaces then 
summing each element matrix during the assembly process, resulting in 




;=i 


C = £c, 


1=1 


k = Y j k, 


1=1 


dFi dFi dFi , 




d8 


The resultant matrix equation is of the form 

M(k)x e + C(k)x e + K(k)x e =F + 

r 

where x e = 


dF(k) dF(k) dF c 

-ir s+ ^r* + ds 8 


(130) 

(131) 

(132) 

(133) 

(134) 


is the nodal displacement vector. 


01 Wl w l Vi Vj ••• w„ 0, I+ 1 w„+iw n+1 V„+l v n+l 
M is the mass matrix, C is the damping matrix, K is the stiffness, and F is the force vector, all as functions of the 
reduced frequency parameter k. 
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B. Structural Damping 


It is important to note that the aerodynamic damping matrix can be augmented by a structural damping matrix. The 
structural damping matrix can be obtained by transforming the generalized coordinates into modal coordinates via the 
eigenvalue analysis as follows: 

Consider the zero-speed structural dynamic equations 


x e +M s l C s x e +M s 1 K s x e =M S 1 


dF dF dF _ \ 


(135) 


where M s is the structural mass matrix, C s is the structural damping matrix, K s is the structural stiffness matrix at zero 
speed (V„ = 0), and F is the force vector. Assuming that the eigenvalues of the matrix M s 1 K s are positive real and 
distinct, the matrix M s 1 K s may be simplified using the similarity transformation 

m; ] k s = x s q.]x;' (136) 


where X s is the eigenvector matrix and Q. s = diag (ft)] .Oh. ... . C0„) is a diagonal matrix whose entries are the frequencies 
of the structural dynamic modes. 

The structural dynamic matrix equation may therefore be expressed as 


Xe T Af y Cs-Xe T 2C s £l s X s x e — M s 


dF dF dF \ 
+ ^ S + ^J S+ d5 5 ) 


Multiplying through by X 1 and letting <p = X l x e be the modal coordinates gives the transformed structural 
dynamics equation 

dF dF dF \ 

F+ ^s S+ ^s S+ d8 8 ) (137) 


< p +X~ 1 M~ l C s X s q)+Q. 2 s (p = X~ x M~ l 


which can also be expressed in modal coordinates as 


9 dfi dfi . dfi 

< Pi + 2^j0)j(pi + COj (Pi = fi+ d 


(138) 


where Q is the viscous damping ratio of the i-th mode and is typically a parameter obtained by ground vibration testing 
and similar methods. 

If £ = d i ag ( £ | , £2 gives the diagonal viscous damping ratio matrix, then the structural damping matrix is 

computed as 

C s = 2 M S X S C Q. S X~ 1 (139) 

The total damping matrix is then given by the linear superposition of the structural and aerodynamic damping 
matrices 

C = C s +C a (140) 

where C a is the aerodynamic damping matrix computed from the previous section. 

The system of equations is then translated into a state space form 


Xe 


0 / 
—M~ X K —M~ X C 


X e 

Xe 



(141) 


with I representing the identity matrix, whereupon the eigenvalues of the matrix equation yield the vibrational fre- 
quencies and mode shapes of the aeroelastic system. The flutter boundary is defined to be the airspeed at which the 
real parts of the eigenvalues of the systems become zero. 


C. Time Integration Methods 

The state space equations may be used to perform time integration of the nodal displacement values given an initial 
deflection profile. For instance, the first-order, explicit (forward) Euler scheme of integration at time step i + 1 is given 

by 



0 

-M- l K 
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-M~ l C 


X e 
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AT 1 



dF < 

ds S 


dF_ 
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S + 




( 142 ) 
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where 


2 


(143) 


At < 
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eig | 


0 I 

- M-'K -M~ l C 


must be satisfied in order to maintain numerical stability of the solution. 

This integration scheme can become prohibitively expensive to enforce in problems with large frequency magni- 
tudes. One approach is to apply modal truncation method so that the time increment may be artificially increased by 
truncating the number of modes obtained from the analysis. This is due to the fact that it is generally difficult to excite 
higher-frequency modes because the energy in a structure is proportional to jKx 2 , which is in turn proportional to 
the squared angular frequency 0 )} r However, there are alternative time integration methods which can be employed 
to relax such restrictions or increase solution accuracy. One such method is the implicit backward Euler scheme of 
integration which is similarly expressed as 
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from which it follows that 


X e 

= { I- At 

Xe 

m V L 


0 

-M- l K 


I 

- M~ l C 


x e 

x e 


T At 

- i 



(144) 


(145) 


This method is fairly stable and permits larger time steps to be used. It is important to note, however, that the 
explicit scheme is time-accurate, as compared to the implicit scheme which effectively damps out higher-frequency 
oscillations. Due to the first-order nature of both the Euler methods, higher-accuracy methods may be desired. One 
popular solution method, Newmark integration, 14 bypasses the state space form and instead algebraically updates the 
nodal degrees of freedom of the aeroelastic finite element equations: 


Mx + Cx + Kx = F 


(146) 


Given specified initial values of x, x, and x (assumed to be zero in this study), one can freely select a time step A t 
and the parameters and /j >i (5 + |)T Then, the integration coefficients are defined as 
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Then, for each time step i, the integration is performed as follows: 


Ff =F{ + M (aoXi + ci2Xi + a^Xi) + C (aw + ci4Xi + a^Xi) (156) 

x i+ 1 =K~ l Fi (157) 

x i+ i =a 0 (x i+ i - x t ) - a 2 Xi - a 3 x,- (158) 

X;+i =Xi + a&Ci + a 2 Xi+\ (159) 


It can be shown that for 8 = j and P = \, the Newmark time integration method is second-order accurate and 
unconditionally stable. However, due to the numerous integration parameters and multi-step nature of the method, the 
Newmark scheme can be computationally expensive. 


V. Flight Dynamic Coupling 

Consider the rigid aircraft flight dynamics in the stability axes described by 

m a ( u + qw — rv) = ( Cl sin a— Co cos a) qooS + T — m a g sin 0 (160) 

m (v + ru — pw) =CyqooS + mg cos 0 sin 0 (161) 

m (w + pv — qu) = {—Cl cos OC — Cd sin a) qooS + mg cos 0 cos 0 (162) 

/y xP f vy q A.; f ^x:Pq 7 vv p/' (/)- Iyy) qf Iy z {l~~ q~ ) =CiqooSb (163) 

^xyP T lyyq A.f T lyzPq (yy'T f ( 7vx C ) PX T I_xz (p~ x ~ ) =C m q !X ,Sc -f TZe (164) 

^xzP lyztf CT ^x:PX f / vv Ixx) PQ ~\~^xy (<7 7* ) C n qoaSb (165) 

0 + g sin 0 tan 0 + rcos 0 tan 0 (166) 

0 =qcos 0 — rsin0 (167) 

=< 7 sin 0 sec 0 + rcos 0 sec 0 (168) 


where is the aircraft mass, S is the aircraft reference wing area, l xx , l yy , I zz , Ixy, I xz , I yz are the aircraft principal 
moments of inertia, c is the mean aerodynamic chord, b is the wing span, is the offset of the thrust line below the 
aircraft CG, and (0,0, y/) are the aircraft Euler angles. 

Note that the aerodynamic coefficients Cl, Cd, C y , Ci, C m , and C„ are influenced by the aeroelastic deflections of 
the aircraft wings. So, the equations of motion of rigid aircraft are coupled with the aeroelastic equations. 

The wing contribution to the aircraft unsteady lift coefficient is evaluated as 

a/- 1 n\ 1 ( tz n \ XXacX G {k) KOC mc C 2 a j /i£n\ 

A C L {k) = -J [c La a ac F {k)+c La — + c Ls 8 J cos A cdx (169) 

This can be also written as 
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Assuming a parabolic drag polar, the aircraft unsteady drag coefficient is evaluated from the unsteady lift coeffi- 
cient as 

C 2 L (k) 

" M 


C D (fc) - C 0o + <5 ( , + C D 5; + C 0 ^ t> r + C Di2 t> ; 2 


(181) 


kARe ' 1 1 

where Cn , C/) 2 , Cd s , and C/j 2 are the drag derivatives due to the elevator and rudder deflections. 

£ Og f Oy 

Neglecting the drag contribution, the wing contribution to the aircraft unsteady pitching moment coefficient is 
evaluated as 
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Due to symmetry, the partial derivative contributions to the aircraft unsteady lift and pitching moment do not 
involve aircraft lateral-directional states /3 , p, and r since these terms cancel out in the integration over both the left 
and right wings. 

The aircraft unsteady rolling moment coefficient is evaluated as 


Cl ( k ) — ££ J ^ (y acC La a acF ( k ) + y ac CL ( 
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The aircraft unsteady yawing moment coefficient is evaluated as 


Cn{k) ~ SbL 
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kARe J 
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(184) 


(185) 


Due to symmetry, the partial derivative contributions to the aircraft unsteady rolling and yawing moments do not 
involve aircraft longitudinal states a and q since these terms cancel out in the integration over both the left and right 
wings. 


VI. Vortex-Lattice Aerodynamic Model Coupling 

Vorview is a computational aerodynamic tool that is used for the development of the aeroelastic computational 
capability. 15 Vorview provides a rapid method for estimating aerodynamic force and moment coefficients as well 
as aerodynamic stability and control derivatives for a given aircraft configuration. It is based on the vortex-lattice 
lifting line aerodynamic theory. The vehicle configuration is constructed within Vorview by a series of panels that 
are formed by spanwise and chordwise locations of bound vortices. Vorview computes the vehicle aerodynamics 
in both the longitudinal and lateral directions independently. The longitudinal and lateral aerodynamics are then 
combined to produce overall aerodynamic characteristics of the vehicle at any arbitrary angle of attack and angle of 
sideslip. Due to the inviscid nature of any vortex-lattice method, the drag prediction by Vorview is most reliable 
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for induced drag prediction. For viscous drag due to boundary layer separation or wave drag due to shock-induced 
boundary layer separation, the prediction may be less reliable. Vorview can provide a rapid estimation of aerodynamic 
derivatives including dynamic derivatives due to angular rates. Owing to the computationally efficient vortex-lattice 
method, aerodynamic derivatives can be estimated in Vorview fairly quickly. A flight dynamic model for a given 
vehicle configuration can be easily developed with Vorview that supplies the model with all necessary aerodynamic 
information for the vehicle. Vorview has been validated by both wind tunnel data 7 as well as the NASA Cart3D tool, 16 
which is a high-fidelity inviscid (Euler) CFD analysis code targeted at analyzing aircraft performance in conceptual 
and preliminary aerodynamic design. In general, both Vorview and Cart3D seem to have similar predictive capabilities 
when compressibility is not a factor. 

Figure 1 1 illustrates an aerodynamic model of the GTM in Vorview. 



Fig. 1 1 - Vorview Aircraft Model 

A. Automated Vehicle Geometry Modeling Tool 

To enable a coupled aeroelastic solution, the aircraft deformed geometry must be generated at each iteration. An 
automated vehicle geometry modeling tool has been developed in MATLAB to update the aircraft deformed geometry. 
The vehicle geometry modeler directly outputs a geometry input file that can be read by Vorview during a solution 
cycle. 


Zg Zy 

A A 



■>x B 




Fig. 12 - GTM Coordinate Systems 

With reference to Fig. 12, the coordinate reference frame (xb^b^zb) defines the Body Station (BS), the Body Butt 
Line (BBL), and the Body Water Line (BWL) of the aircraft, respectively. The coordinate reference frame (xy.yy,z.y) 
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is the translated coordinate system attached to the nose of the aircraft. The stability reference frame (x,y,z) is attached 
to the CG such that x = xy — xy , y = y v — yv , and z = zv — zy , where (xy ,yv,zv) is the coordinate of the CG in the 
(. xy,yy,zv ) reference frame. 7 

The vehicle geometry modeler has access to the outer mold line of the jig-shape (undeformed) aircraft geometry. 
The coordinate reference frame (xy,yy,zv) defines the coordinate system used in the vehicle geometry model. 

The aeroelastic deflections in bending and torsion result in a displacement Ap and rotation angle <j> where 

0 = ©d i - W x d 2 + V v d 3 (186) 

Ar = — (W sin W x + V sin V x ) di + V cos V x d 2 + W cos W x d 3 (187) 

The coordinate reference frame (. x,y,z ) is related to the coordinate reference frame {xy,yy,zy) by the following 
relationship: 
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sinAsinT 

cos A sin r 
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(188) 


where (v i , V 2 , v 3 ) are the unit vectors for the Vorview coordinate reference frame (xy ,yv,zv). 

Thus, the aeroelastic deflections result in a wing twist expressed as an incremental rotation vector (A0*, A0y, A0 Z ) 
and a displacement vector (Axy,Ayy,Azy) 


A 0c =®sinAcosr — W t cosA — y t sinAsinr (189) 

A0 V = — ©cosAcosT — iy t sinA + y v cosAsinr (190) 

A0 ; =©sinr + y t cosr (191) 

A xy = — (W sin W x + V sin V x ) sin A cos T + V cos V x cos A — W cos W x sin A sin T (1 92) 

A yy =(W sinW t + y sin V x ) cos A cos T + V cosy v sinA + W cos W x cos A sin F (193) 

Azy = — (W sin W x + V siny v ) sinT + W cosiy T cosr (194) 


A coordinate transformation to account for wing aeroelastic deflections is performed by rotating a wing section 
about its elastic axis by the incremental rotation vector (A<j> x , A<j) y , A<j> z ) and then translating the resultant coordinates 
by the displacement vector (Axy, Ayy, Azy). 

To perform the coupled aeroelastic computation, the static aeroelastic model is coupled with Vorview through the 
automated vehicle geometry modeler. Aerodynamic force and moment coefficients as computed from Vorview are used 
as inputs to the static aeroelastic model. The computed aeroelastic deflections are then used to generate the aircraft 
deformed geometry by the automated vehicle geometry modeler. The aerodynamic solution is then recomputed with 
the aircraft deformed geometry in Vorview. This process is iterated until the solution is converged when errors in the 
computed aeroelastic deflections are within a specified tolerance. A flow chart for the coupled aeroelastic computation 
is shown in Fig. 13. 

The static aeroelastic solution provides aerodynamic information for the deformed aircraft under a trimmed flight 
condition. The dynamic aeroelastic analysis is conducted to compute the unsteady contributions to wing aerodynamics. 
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Fig. 13- Coupled Aeroelastic Vortex Lattice Computation Flow Chart 


VII. Simulations 


A coupled aeroelastic-longitudinal flight dynamic model is built by coupling the wing dynamic aeroelastic model 
to a linearized model of the aircraft rigid-body longitudinal flight dynamics. The wing dynamic aeroelastic model is 
represented by a second-order system described by Eq. (134) which is assembled using the finite-element method. 
In this implementation, the wing dynamic aeroelastic model is implemented using the cantilever boundary conditions 
where W (0) , VT v (0), V (0), V't(O) , and 0(0) are all set to zero. Strictly speaking, the wing symmetric modes should be 
coupled with the aircraft longitudinal flight dynamics where the wing displacement at the wing root matches the aircraft 
displacement. By removing the wing displacement at the wing root through the cantilever boundary conditions, the 
coupled aeroelastic-longitudinal flight dynamic model is a reasonable approximation of the coupled wing symmetric 
modes with the aircraft rigid-body modes. 

Let the aircraft rigid-body flight dynamics be represented by a first-order system given by 


where x r is the aircraft rigid-body state vector, u 


M r x r = Sx r + C u u 
r -i T 


(195) 



8 


is a control vector of the control surface deflections 


comprising the elevator 8 e and the symmetric VCCTEF deflection 8 , and C„ 


Cg e C , 5 


is the partitioned control 


sensitivity matrix. 

Note that the design concept of the VCCTEF utilizes only the third camber segments of each spanwise flap section 
for roll and pitch control due to the faster response of these control surfaces provided by the electric drive motors. 

r -\ T 


Thus, 8 = 


<5i <>2 83 ... S l5 


represents the deflections of the third camber segments of the 15 spanwise 


flap sections of the VCCTEF. 

Due to the coupling between the aircraft unsteady aerodynamic coefficients with the aeroelastic deformation of the 
wing, the coupling matrices represented by H\ (k). H 2 (k). and Hj, ( k ) are introduced, where k is the reduced frequency 
parameter. These matrices are calculated using the unsteady aircraft aerodynamic coefficients evaluated using the 
equations developed in the previous section. Thus, the first-order system representing the aircraft rigid-body dynamics 
coupled to the wing aeroelastic states is given by 


M, x r = Sx r + H 1 ( k)x e + H 2 (k)x e + Hj, ( k)x e +C u u 


(196) 
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The coupled equations (196) and (134) can be combined together to form a coupled first-order state space model 
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Note that the coupled state space model represents a reduced-frequency-dependent state space model. This model 
is generally valid for a known value of the reduced frequency parameter k and can be used to approximate a flight dy- 
namic model if the dominant wing aeroelastic frequency is known. Another method is to remove the reduced frequency 
dependency from the state space model by approximating the Theodorsen’s function C{k) by various methods such as 
the Roger method of rational fraction approximation 17 or the R.T. Jones method. 18 These approximation methods can 
be advantageous in that the state space model is more broadly applicable to a wider range of frequencies at a cost of 
increasing the order of the state space model by introducing aerodynamic lag states resulting from the approximation 
of the Theodorsen’s function. This approach will be considered in the future work. 


A. Linearized Aircraft Rigid-Body Longitudinal Flight Dynamic Model 


A 4-state longitudinal flight dynamic model of the GTM is implemented. The aircraft rigid-body states are given by 
lT 


X r = 


f Aa 


q A 6 


where 


AV 


v is a the normalized perturbation of the aircraft forward airspeed, Aa is the 
perturbation of the aircraft angle of attack, q is the pitch rate, and Ad is the perturbation in the aircraft pitch angle. For 
the flight condition of Mach 0.8 at 35,000 ft, the matrices for the aircraft rigid-body longitudinal flight dynamic model 
are given by 



11.1138 

0 

0 

0 


0 

11.1757 

0 

0 


0 

0.1310 

0.7841 

0 


0 

0 

0 

1 


-0.0558 

-0.4364 

-0.7480 

-0.4595 

-1.7284 

-6.3068 

10.9544 

-0.0306 

-0.0074 

-1.7648 

-0.3370 

0 

0 

0 

1 


0 


which represent a linearization about the rigid-body trim point a = 3.8142°, V = 778.2063 ft/sec, 6 = —3.8142°, 
T = 5617 lbs, and 8 e = -6.1497°. 

The eigenvalues of the aircraft rigid-body longitudinal flight dynamics are calculated to be 

X sp =-0.5779 ±1.449 If 
X p = - 0.0042 ± 0.0763; 

The eigenvalues X sp correspond to the high frequency, highly damped short-period mode of the aircraft dynamics. 
The eigenvalues X p represent the low frequency, lightly damped phugoid mode of the aircraft dynamics involving the 
airspeed, pitch angle, and altitude. 

The natural frequencies and the damping ratios for the short -period and the phugoid modes are calculated to be 


a>n sp =1.5601 rad/sec 
C sp =0.3704 
a )„ p =0.0764 rad/sec 
Cp =0.0545 
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B. Flutter Analysis 

Two different wing models will be analyzed for the baseline wing stiffness of the GTM wings and for the reduced 
stiffness of the highly flexible wings of the ESAC. The baseline structural rigidities El and GJ are estimated for the 
conventional stiff wing structures for the GTM. For the ESAC, the wing structural rigidities El and GJ are purposely 
reduced by a factor of two to model highly flexible wing structures. The increased flexibility enables the wing shaping 
control actuation by the VCCTEF system for drag reduction. 

In addition to the wing dry mass, the fuel mass is also accounted for. The fuel is stored in the center tank and 
wing main tanks. The center tank holds 20,000 lbs of fuel. Each of the main tanks holds about 15,000 lbs of fuel. 
The center tank is used first until it is empty. Then, the fuel is drawn equally from the wing main tanks. The fuel 
mass is modeled as the combined wing mass density. As the structural rigidities are reduced, the wing dry mass also 
decreases. Assuming that the wing box structure is modeled as a thin-walled structure, then the mass change is related 
to the change in the wing structural rigidity El can be modeled. 

For the flutter analysis, the structural dynamic modes for the cantilever, symmetric, and anti-symmetric boundary 
conditions are computed with 80% fuel loading and no structural damping for conservatism. A trim thrust value is 
used in the flutter analysis based on the linearization of the aircraft rigid-body longitudinal flight dynamic model to 
account for aero-propulsive-elastic effects. The structural dynamic cantilever mode shapes of the stiff GTM wings 
are shown in Fig. 14 and their associated natural frequencies for the stiff GTM wings and flexible ESAC wings are 
summarized in Table 1 . 


Mode 1 B 


Mode 2B 








Fig. 14 - Structural Dynamic Cantilever Mode Shapes of Stiff GTM Wings 


Mode 

Natural Frequency, Hz (GTM Wings) 

Natural Frequency, Hz (ESAC Wings) 

IB 

1.4934 

1.1252 

2B 

4.0522 

2.9620 

1T/2B 

5.0930 

3.7895 

3B 

9.5070 

7.1341 

2T/3B 

10.9926 

8.4271 

2T/4B 

17.4842 

13.3635 


Table 1 - Structural Dynamic Natural Frequencies of Cantilever Modes with 80% Fuel Loading 
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The flutter analysis using a linear aeroelastic model is conducted for the cantilevered wing with and without 
coupling to the aircraft rigid-body longitudinal flight dynamic model. The flutter analysis of the coupled aeroelastic - 
longitudinal flight dynamic model allows for solving for the flutter frequency that will be used to calculate the reduced 
frequency parameter k needed to generate the coupled state space model. A cruise condition at 35,000 ft altitude is 
examined. 

1. Stiff GTM Wings 

The frequency and damping ratios of the stiff GTM wings with the baseline stiffness were computed by sweeping over 
a Mach number range. Figure 15 is a plot of the aeroelastic frequencies and damping ratios for the uncoupled wing 
cantilever modes, while Table 2 summarizes the critical flutter mach numbers and flutter frequencies for the first two 
flutter modes. The critical flutter mode is observed to be due to the second bending mode (2B) at Mach 1.3792 with a 
flutter frequency of 2.4792 Hz. 

Stiff Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Stiff Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 




Fig. 15 - Flutter Results for Uncoupled Cantilever Modes of Stiff GTM Wings at 35,000 ft with 80% Fuel Loading 


Stiff Wing Coupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Stiff Wing Coupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 



Fig. 16 - Flutter Results for Coupled Cantilever Modes of Stiff GTM Wings at 35,000 ft with 80% Fuel Loading 


Mode 

Flutter Mach 

Flutter Frequency, Hz 

2B 

1.3792 

2.4792 

3B 

1.6729 

6.4578 


Table 2 - First Two Flutter Speeds of Uncoupled Cantilever Modes of Stiff GTM Wings at 35,000 ft with 80% Fuel 

Loading 
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The flutter analysis is repeated for the coupled aeroelastic -longitudinal flight dynamic model. The critical flutter 
mode is observed to be due to the third bending mode (3B) at Mach 1.5490 with a flutter frequency of 9.6130 Hz. 
Note that the flutter characteristics are significantly changed as a result of the coupling with the aircraft rigid-body 
longitudinal flight dynamics. 

2. Flexible ESAC Wings 

The stiffness of the flexible ESAC wings is reduced by half from the baseline stiffness of the GTM wings. The 
flexibility will allow the VCCTEF to be more effective in wing shaping control for drag reduction. As a result of 
the reduced stiffness, the flutter boundary of the flexible ESAC wing will decrease from that of the stiff GTM wings. 
Without coupling to the aircraft rigid-body longitudinal flight dynamics, the flutter characteristics of the flexible ESAC 
wings are shown in Fig. 17 and Table 3. The critical flutter mode is observed to be due to the second bending mode 
(2B) at Mach 1.01 15 with a flutter frequency of 3.2367 Hz. 


Flexible Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Flexible Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 




Fig. 17 - Flutter Results for Uncoupled Cantilever Modes of Flexible ESAC Wings at 35,000 ft with 80% Fuel 

Loading 

Flexible Wing Coupled Cantilever Modes @ 35,000 ft w / 80% Fuel Loading Flexible Wing Coupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 




Fig. 18 - Flutter Results for Coupled Cantilever Modes of Flexible ESAC Wings at 35,000 ft with 80% Fuel Loading 


Mode 

Flutter Mach 

Flutter Frequency, Hz 

2B 

1.0115 

3.2367 

1T/2B 

1.1888 

5.0131 


Table 3 - First Two Flutter Speeds Uncoupled Cantilever Modes of Flexible ESAC Wings at 35,000 ft with 80% Fuel 

Loading 
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The flutter characteristics of the coupled cantilever modes of the flexible wings are shown in Fig. 18. The critical 
flutter mode is observed to be due to the 2T/4B mode which occurs at Mach 1.3105 with a flutter frequency of 6.2244 

Hz. 

C. Modal Analysis 

The flutter analysis can be used to establish the most dominant aeroelastic mode at a given flight condition. This can 
be determined by examining the aeroelastic mode with the lowest damping. In all cases, the second bending mode is 
the most lightly damped mode. The frequency of this mode is then used to determine the reduced frequency parameter 
k for the coupled aeroelastic-longitudinal flight dynamic model. 

1. Stiff GTM Wings 

A modal analysis of the coupled aeroelastic-longitudinal flight dynamics is performed to determine the effect the 
aeroelastic coupling on the aircraft rigid-body modes. The eigenvalues of the aircraft rigid-body longitudinal flight 
dynamics with the stiff GTM wings are calculated to be 

X sp =-0.6012 ±1. 3949/ 

X p =-0.0051 ± 0.0771/ 

The natural frequencies and the damping ratios for the short -period and the phugoid modes are calculated to be 

(On sp =1.5 189 rad/sec 
=0.3958 

(O np =0.0773 rad/sec 
C p =0.0655 

The frequencies and damping ratios for the coupled cantilever modes of the stiff GTM wings are shown in Table 4. 
The root locus of the poles for the phugoid, short-period, and the first three cantilever modes of the stiff GTM wings 
are plotted in Fig. 19. 


Mode 

Eigenvalue 

Frequency (rad/s) 

Damping Factor 

IB 

-0.7755 ±10.4433/ 

1.6667 

0.0741 

2B 

-0.3637 ±25.8867/ 

4.1204 

0.0140 

1T/2B 

-0.6130 ±32.2080/ 

5.1270 

0.0190 

3B 

-0.7252 ±59.5928/ 

9.4852 

0.0122 

2T/3B 

-1.2371 ±66.9549/ 

10.6580 

0.0185 

2T/4B 

-1.5012 ±108.3448/ 

17.2453 

0.0139 


Table 4 - Natural Frequencies and Damping Ratios of Coupled Cantilever Modes of Stiff GTM Wings at 35,000 ft 

with 80% Fuel Loading 


32 of 38 


American Institute of Aeronautics and Astronautics 




Root Locus of Stiff Wing Cantilever Modes @ 35,000 ft w / 80% Fuel Loading 
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Fig. 19 - Root Locus of Coupled Aeroelastic -Longitudinal Flight Dynamics with Stiff GTM Wings at 35,000 ft with 

80% Fuel Loading 


2. Flexible ESAC Wings 

The eigenvalues of the aircraft rigid-body longitudinal flight dynamics with the flexible ESAC wings are calculated to 
be 


X sp =- 0.6570 ±1.2712/ 

X p =- 0.0063 ± 0.0782/ (198) 

The natural frequencies and the damping ratios for the short -period and the phugoid modes are calculated to be 

(0„ sp =1.43 10 rad/sec 
C s P =0.4591 
( 0„ p =0.0785 rad/sec 
c p =0.0807 

Note that both the rigid-body modes are stable in spite of the aeroelastic coupling. This is due to the significant 
frequency separation between the rigid-body modes and the wing cantilever modes. 

The frequencies and damping ratios for the coupled cantilever modes of the stiff GTM wings are shown in Table 5. 
The root locus of poles for the phugoid, short-period, and the first three cantilever modes of the flexible ESAC wings 
are plotted in Fig. 20. 


Mode 

Eigenvalue 

Frequency (rad/s) 

Damping Factor 

IB 

-1.0929 ±9. 3646/ 

9.8499 

0.1159 

2B 

-0.2881 ±19.6302/ 

19.6323 

0.0147 

1T/2B 

-1.0147 ±25.3135/ 

25.3339 

0.0401 

3B 

-0.8460 ±46.4797/ 

46.4874 

0.0182 

2T/3B 

-5.1495 ±46.6200/ 

46.9035 

0.1098 

2T/4B 

-6.7807-81.6900/ 

81.9709 

0.0827 


Table 4 - Natural Frequencies and Damping Ratios of Coupled Cantilever Modes of Flexible ESAC Wings at 35,000 

ft with 80% Fuel Loading 
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Root Locus of Flexible Wing Cantilever Modes @ 35,000 ft w / 80% Fuel Loading 
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Fig. 20 - Root Locus of Coupled Aeroelastic-Longitudinal Flight Dynamics with Flexible ESAC Wings at 35,000 ft 

with 80% Fuel Loading 


D. Dynamic Response Simulations 

1. Elevator Input 

The dynamic response for aircraft with the stiff GTM wings with the baseline stiffness is simulated for separate inputs 
of the elevator and VCCTEF over a time span of 40 sec. Figure 21 shows a doublet input of the elevator. 


1.5 | 1 ! 1 . ! ! ! 

1 y — : : : : : - 

0.5 - -p - 

0 - i— i— 

-o.5 - ; - ; I — ; — ■ ■■; - 

-1 - f i •: - 

-1.5 1 1 1 1 1 1 1 1 

0 5 10 15 20 25 30 35 40 

t, sec 


Fig. 21 - Elevator Doublet 

The dynamic responses of the aircraft with rigid-body and coupled aeroelastic-longitudinal flight dynamics with 
the stiff GTM wings and flexible ESAC wings are shown in Fig. 22. The dynamic response due to the stiff GTM 
wings matches that of the rigid-body longitudinal flight dynamics very well. The dynamic response due to the flexible 
ESAC wings also matches the rigid-body dynamic response reasonably well, but some small differences are noted. 
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Fig. 22 - Aircraft Longitudinal States 

The dynamic responses of the aeroelastic deflections at the wing tip, denoted as W np and © r , p are shown in Fig. 
23. The flexible ESAC wings experience much greater aeroelastic deflections as expected since the stiffness of the 
ESAC wings is half that of the stiff GTM wings. It is somewhat surprising that, with the significant wing aeroelastic 
deflections, the aircraft longitudinal states do not seem to be much affected. 
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Fig. 23 - Aeroelastic Deflections at Wing Tip 
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2. VCCTEF Input 


The VCCTEF can be used for either roll or pitch control by symmetric or anti-symmetric deflections of the individual 
spanwise flap segments. The first inboard flap segment is designated as a high-lift flap and therefore is not used for 
roll or pitch control. Because of the continuous trailing edge surfaces, the flap segment adjacent to the inboard flap 
can only be deflected by a relative amount as permitted by the transition material. 

For the dynamic response simulation of the VCCTEF, the outboard flap, designated flap 15, is commanded by the 
same doublet as the elevator. The commands for flaps 1 to 15 vary linearly from zero to the full doublet. 

The dynamic responses of the aircraft with rigid-body and coupled aeroelastic-longitudinal flight dynamics with 
the stiff GTM wings and flexible ESAC wings are shown in Fig. 24. These dynamic responses are significantly 
different from one another. In particular, the dynamic response due to the flexible ESAC wings exhibits control 
reversals for a, q , and 0. The control reversals are due to a large nose-down twist caused by the flexible ESAC wings. 






Fig. 24 - Aircraft Longitudinal States 

The dynamic responses of the aeroelastic deflections at the wing tip, denoted as W t i p and © r , p are shown in Fig. 
25. The flexible ESAC wings experience much greater aeroelastic deflections and dynamic transients than the stiff 
GTM wings. In particular, the twist of the flexible ESAC wings is quite substantial, more than 2 degrees at some time 
instances. This could explain the sign reversal in a , q , and 9. 
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Fig. 25 - Aeroelastic Deflections at Wing Tip 

VIII. Conclusions 

This paper presents a coupled aeroelastic-longitudinal flight dynamic model of a flexible wing aircraft. This aircraft 
concept called Elastically Shaped Aircraft Concept (ESAC). The aircraft concept addresses the drag reduction goal in 
commercial aviation through an elastic wing shaping control approach for aircraft with highly flexible wing structures. 
The multi-disciplinary nature of flight physics is appreciated with the recognition of the potential adverse effects of 
aeroelastic wing shape deflections on aerodynamic performance. By aeroelastically tailoring the wing shape with 
active control, a significant drag reduction benefit could be realized. To attain the potential of the elastic wing shaping 
control concept, a new type of aerodynamic control effector is introduced and is referred to as a Variable Camber 
Continuous Trailing Edge Flap (VCCTEF). 

A coupled aeroelastic -flight dynamic modeling approach has been developed for this aircraft concept. The flight 
dynamic model is coupled with the aeroelastic states from a finite-element model of the flexible wing via the aeroe- 
lastic contributions to the aerodynamic coefficients. A coupled aeroelastic-longitudinal flight dynamic model has 
been developed for both the stiff wing GTM and flexible wing ESAC. Initial simulations show that the short-period 
and phugoid modes remain stable for the flexible wing ESAC. An open-loop response simulation is conducted to 
demonstrate the coupled dynamic response. The wing flexibility results in a significant deflection as compared to a 
conventional stiff wing transport aircraft which could cause some issues of control reversal. Future work will further 
develop the coupled aeroelastic lateral-directional dynamic model and ultimately a fully coupled 6-dof flight dynamic 
model. 
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